{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage import data\n",
    "camera = data.camera()\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<span style=\"font-size: 150%;\">**Exercise:** Write an algorithm from scratch that will\n",
    "take an input image as an ndarray, and rotate it by $\\theta$ degrees.</span>\n",
    "\n",
    "If you feel creative, you can also write code to magnify (zoom) the image.\n",
    "<p></p>\n",
    "You may need: http://en.wikipedia.org/wiki/Polar_coordinate_system\n",
    "<p></p>\n",
    "A (bad) solution is given below--but try it yourself before looking!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def interp(x, y, vals):\n",
    "    \"\"\"Bilinear interpolation.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    x, y : float\n",
    "        Position of required value inside a 1x1 pixel.\n",
    "    vals : list of 4 floats\n",
    "        Values on the four corners of the pixels, in order\n",
    "        top left, top right, bottom left, bottom right.\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    v : float\n",
    "        Interpolated value at position ``(x, y)``.\n",
    "    \n",
    "    \"\"\"\n",
    "    top_left, top_right, bottom_left, bottom_right = vals\n",
    "    \n",
    "    top_middle = (1 - x) * top_left + x * top_right\n",
    "    bottom_middle = (1 - x) * bottom_left + x * bottom_right\n",
    "    \n",
    "    return (1 - y) * top_middle + y * bottom_middle\n",
    "\n",
    "\n",
    "def rotate(image, theta):\n",
    "    theta = np.deg2rad(theta)\n",
    "    \n",
    "    height, width = image.shape[:2]\n",
    "    out = np.zeros_like(image)\n",
    "    \n",
    "    centre_x, centre_y = width / 2., height / 2.\n",
    "    \n",
    "    for x in range(width):\n",
    "        for y in range(height):\n",
    "            \n",
    "            x_c = x - centre_x\n",
    "            y_c = y - centre_y\n",
    "            \n",
    "            # Determine polar coordinate of pixel\n",
    "            radius = np.sqrt(x_c**2 + y_c**2)\n",
    "            angle = np.arctan2(y_c, x_c)\n",
    "            \n",
    "            new_angle = angle - theta\n",
    "            \n",
    "            old_x = radius * np.cos(new_angle)\n",
    "            old_y = radius * np.sin(new_angle)\n",
    "            \n",
    "            old_x = old_x + centre_x\n",
    "            old_y = old_y + centre_y\n",
    "            \n",
    "            if (old_x >= width - 1) or (old_x < 1) or\\\n",
    "               (old_y >= height - 1) or (old_y < 1):\n",
    "                    continue\n",
    "            else:\n",
    "                xx = int(np.floor(old_x))\n",
    "                yy = int(np.floor(old_y))\n",
    "                \n",
    "                out[y, x] = interp(old_x - xx, old_y - yy,\n",
    "                                   [image[yy, xx], image[yy, xx + 1],\n",
    "                                    image[yy + 1, xx], image[yy + 1, xx + 1]])\n",
    "    \n",
    "    return out\n",
    "\n",
    "rotated = rotate(camera, 40)\n",
    "    \n",
    "plt.imshow(rotated, cmap='gray', interpolation='nearest');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Advanced challenge: rectifying an image\n",
    "\n",
    "<img src=\"../../images/chapel_floor.png\" style=\"float: left;\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rectify the tiles above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from skimage import transform\n",
    "\n",
    "\n",
    "from skimage.transform import estimate_transform\n",
    "\n",
    "source = np.array([(129, 72),\n",
    "                   (302, 76),\n",
    "                   (90, 185),\n",
    "                   (326, 193)])\n",
    "\n",
    "target = np.array([[0, 0],\n",
    "                   [400, 0],\n",
    "                   [0, 400],\n",
    "                   [400, 400]])\n",
    "\n",
    "tf = estimate_transform('projective', source, target)\n",
    "H = tf.params   # in older versions of skimage, this should be\n",
    "                # H = tf._matrix\n",
    "\n",
    "print(H)\n",
    "\n",
    "# H = np.array([[  3.04026872e+00,   1.04929628e+00,  -4.67743998e+02],\n",
    "#               [ -1.44134582e-01,   6.23382067e+00,  -4.30241727e+02],\n",
    "#               [  2.63620673e-05,   4.17694527e-03,   1.00000000e+00]])\n",
    "\n",
    "def rectify(xy):\n",
    "    x = xy[:, 0]\n",
    "    y = xy[:, 1]\n",
    "\n",
    "    # You must fill in your code here.\n",
    "    #\n",
    "    # Handy functions are:\n",
    "    #\n",
    "    # - np.dot (matrix multiplication)\n",
    "    # - np.ones_like (make an array of ones the same shape as another array)\n",
    "    # - np.column_stack\n",
    "    # - A.T -- type .T after a matrix to transpose it\n",
    "    # - x.reshape -- reshapes the array x\n",
    "\n",
    "    # We need to provide the backward mapping\n",
    "    HH = np.linalg.inv(H)\n",
    "\n",
    "    homogeneous_coordinates = np.column_stack([x, y, np.ones_like(x)])\n",
    "    xyz = np.dot(HH, homogeneous_coordinates.T)\n",
    "\n",
    "    # We want one coordinate per row\n",
    "    xyz = xyz.T\n",
    "\n",
    "    # Turn z into a column vector\n",
    "    z = xyz[:, 2]\n",
    "    z = z.reshape([len(z), 1])\n",
    "\n",
    "    xyz = xyz / z\n",
    "\n",
    "    return xyz[:, :2]\n",
    "\n",
    "image = plt.imread('../../images/chapel_floor.png')\n",
    "out = transform.warp(image, rectify, output_shape=(400, 400))\n",
    "\n",
    "f, (ax0, ax1) = plt.subplots(1, 2, figsize=(8, 4))\n",
    "ax0.imshow(image)\n",
    "ax1.imshow(out)\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1+"
  },
  "widgets": {
   "state": {},
   "version": "1.1.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
